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Statement  of  Progress: 

Theoretically  we  found  that  when  a non- ideal  blast  wave  and  exothermic 
reaction  kinetics  are  modeled  realistically,  both  the  total  source  energy 
and  the  source  energy  density  are  important  to  the  direct  initiation  process. 

At  high  energy  densities  the  total  energy  is  limiting,  and  as  the  energy 
density  decreases  the  total  energy  required  increases  until  at  some  critical 
low  energy  density  the  energy  density  itself  becomes  the  limiting  factor  for 
initiation.  This  confirms  the  experimental  observations  made  by  John  Lee  at 
McGill  University.  We  have  also  learned  how  the  non-ideality  of  a source  re- 
gion influences  the  structure  of  a blast  wave.  Specifically,  the  blast  wave 

from  bursting  spheres,  the  ramp  addition  of  energy,  constant  velocity  centrally 

» 

ignited  spherical  flames,  and  accelerated  spherical  flames  have  been  studied  " 
extensively,  and  this  has  yielded  the  generalizi .tion  that  far  field  positive 
impulse  is  properly  scaled  by  the  source  energy  but  that  for  low  power  density 
sources  far  field  overpressure  falls  below  the  expected  energy  scaled  value. 

Experimentally  the  initiation  kinetics  of  propylene  oxide-nitrogen- 
oxygen  mixtures  was  studied  over  the  temperature  range  of  1000-2500°K.  These 
results  show  the  usual  Arrhenius  dependence  on  temperature  for  the  ignition 
delay. 

This  work  has  led  to  three  M.S.  and  two  Ph.D.  theses,  two  interim  tech- 
nical reports,  four  papers  in  the  open  literature  and  seven  presentations  at 
meetings  and  seminars. 
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ABSTRACT 

This  study  is  both  a theoretical  and  experimental  investigation  of 
the  initiation  of  detonation  by  non-ideal  blast  waves.  The  theoretical 
portions  consisted  of  two  studies:  (1)  an  analytic  study  of  the  flow  fields 
associated  with  heat  addition  during  the  initial  phases  of  heat  addition  to 
a source  region  which  is  chemically  reactive,  and  (2)  a study  of  the  initia- 
tion of  explosions  by  non-ideal  sources  using  a numerical  technique.  The 
analytic  portion  of  this  study  uncovered  many  difficulties  with  the  handling 
of  the  heat  addition  in  a reactive  source  region.  These  are  discussed  ex- 
tensively in  this  report.  The  numerical  calculation  showed  that  when  the 
source  region  was  a bursting  sphere  of  high  energy  density,  the  total  energy 
in  the  sphere  was  the  limiting  factor  for  initiation  and  that  when  the  energy 
density  of  the  sphere  was  low  the  energy  density  was  the  limiting  factor  for 
initiation.  Furthermore,  as  the  energy  density  dropped  the  energy  required 
for  initiation  increased  markedly.  This  agrees  quite  well  in  a qualitative 
sense  with  experimental  observations  from  spark  initiated  detonations.  The 
experimental  program  also  consisted  of  two  phases.  In  one  of  these,  reflected 
initiation  of  propylene  oxide-nitrogen-oxygen  mixtures  was  studied  over  the 
temperature  range  of  about  1000-2500°K.  The  data  was  correlated  to  an  Arrhen- 
ius temperature  dependency  and  yielded  an  activation  energy  of  approximately 
38,000  cal/gm  mole.  The  other  experimental  study  investigated  the  effect  of 
inhibitors  and  promotors  on  the  delayed  time  to  ignition  for  propane-air 
mixtures,  also  using  the  reflected  shock  technique. 
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DIRECT  INITIATION  OF  DETONATION 

I INTRODUCTION 

The  program  of  research  reported  here  was  directed  at  two  Air  Force 
needs,  namely  the  easy  evaluation  of  prospective  fuels  for  fuel-air  ex- 
plosives (FAE)  devices  as  well  as  the  evaluation  of  fuel  tank  vulnerability 
to  ignition  combustion  and  detonation  from  the  penetration  of  hot  multi- 
source  gunfire  or  missile  fragments.  The  common  factor  for  these  two  needs 
is  that  they  depend  upon  the  direct  initiation  of  detonation.  This  phrase 
means  that  process  during  which  an  igniter  causes  very  rapid  transition 
from  the  initial  quiescent  state  of  a reactive  mixture  to  the  final  ongoing 
detonation  in  the  reactive  mixture.  Hence,  if  we  are  to  comprehend  those 
events  which  must  occur  during  the  transition  process,  we  shall  need  to  know 
what  exists  at  the  end  of  the  transition  process.  Thus,  we  first  describe 
the  structure  of  a self-sustaining  propagating  detonation  wave  in  considerable 
detail  and  then  describe  our  current  understanding  of  initiation  behavior. 

This  is  followed  by  sections  devoted  to  details  of  the  work  performed  under 
the  grant. 

II  DETONATION  STRUCTURE 

It  is  well  known  that  a rather  simple,  steady-flow  hydrodynamic 
theory  (the  Chapman-Jouguet  or  CJ  theory)  yields  a reasonable  description  of 
the  gross  properties  of  self-sustaining  detonation  waves. ^ The  waves  in- 
deed appear  to  be  constant  velocity  waves  traveling  at  a high  supersonic 
velocity  as  predicted  by  the  theory  and  they  contain  average  shock  and  CJ 


plane  pressures  which  are  quite  adequately  predicted  by  the  theory.  How- 
ever, the  theory  has  one  rather  severe  limitation  in  that  it  predicts  that 
all  exothermic  systems  should  exhibit  detonation  behavior  and,  therefore, 
that  there  should  be  no  concentration  or  pressure  limits  to  the  propagation 
of  detonation  waves.  Experimentally,  we  know  that  this  result  is  not  true. 
Detonations  in  tubes  do  exhibit  failure  behavior  both  as  the  mixture  is 
diluted  with  inerts  and  as  the  initial  pressure  is  reduced.  Also,  unconfined 
detonations  definitely  exhibit  concentration  limits. 

It  is  now  known  that  the  CJ  concept  only  applies  to  the  gross  overall 
behavior.  If  one  looks  at  the  wave  structure  in  sufficient  detail,  one  al- 
ways finds  that  the  wave  is  actually  propagating  with  a three-dimensional 
non-steady  structure. ^ This  behavior  in  no  way  negates  the  utility  of  the 
CJ  concept  in  describing  the  overall  (or  gross)  behavior  of  the  wave  but 
simply  means  that  we  now  know  that  the  CJ  concept  cannot  be  used  to  look  at 
detailed  structure,  just  as  we  already  knew  that  it  could  not  be  used  to 
discuss  the  detailed  fluid  behavior  near  the  front. 

In  the  past  five  to  ten  years  we  have  learned  a great  deal  about  the 
structure  of  gas  phase  detonation  waves.  We  now  know  that  a one-dimensional 
steady  wave  should  not  exist  because  it  is  hydrodynamical ly  unstable  to  a 
perturbation. v We  also  have  learned  that  the  instability  forms  transverse 
perturbation  waves  traveling  across  the  front  behind  the  main  (or  lead)  shock 
wave  of  the  detonation  in  the  flow  region  which  contains  the  exothermic  chemi- 
cal reactions/3** ^Further,  we  know  that  in  an  "overall  steady"  self-sus- 
taining detonation  these  waves  are  pressure  waves  of  finite  amplitude  and 
propagate  with  an  inherent  spacing  or  structural  size  which  is  related  to  the 
rate  of  the  chemical  reactions  that  are  occurring  in  the  detonation  itself. ^ 
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Since  tho  transverse  waves  are  of  finite  amplitude,  they  are  actually 
reasonably  weak  shock  waves  (Mach  numbers  below  1.4)  which  interact  with 
the  lead  shock  wave  to  produce  Mach  stems  at  their  confluence  with  the  front. 
Furthermore,  because  of  tho  varying  pressure  behind  the  lead  shock  that  is 
caused  by  the  presence  of  these  waves,  the  lead  shock  does  not  have  a con- 
stant velocity.  The  central  element  of  it  typically  shows  a sawtooth-like 
velocity  behavior,  with  the  amplitude  varying  about  ±20  percent  of  the  CJ 
velocity.  This  particular  behavior  occurs  while  each  element  of  the  lead 
shock  wave  is  decaying  at  all  times  with  the  transverse  waves  continually 
overriding  these  decaying  shock  elements.  When  the  transverse  waves  collide, 
they  produce  the  central  shock  element  whose  velocity  jumps  to  a value  well 
above  CJ  immediately  after  the  intersection  of  the  transverse  waves. ^ 

' The  paths  of  the  intersections  of  the  transverse  waves  with  the  lead 
shocks  as  these  shocks  propagate  produce  what  is  called  a detonation  cell. 

A new  cell  starts  with  the  crossing  of  oppositely  directed  transverse  waves 
and  the  cell  ends  at  the  crossing  of  the  next  pair  of  transverse  waves  in  the 
two  trains  of  these  waves. 

Generally,  the  reactions  just  behind  the  lead  shock  are  endothermic  or 
thermally  neutral  before  they  become  exothermic.  This  distance  is  called  the 
induction  distance  and  it  depends  strongly  upon  tho  temperature,  being  shorter 
the  higher  the  temperature  or  the  higher  the  shock  Mach  number.  There  is  an 
upper  limit  to  the  lead  shock  Mach  number,  which  occurs  when  M is  of  the  order 
of  1.5  MCJ,  above  which  the  reactions  do  not  become  exothermic. 

Because  the  induction  distance  is  finite  behind  the  Mach  stem  shock 
created  at  the  start  of  a new  cell,  a hot  spot  is  created  when  the  exothermal 
reactions  occur.  This  hot  spot  generates  a pressure  pulse  which  then 
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overtakes  the  Mach  stem  in  the  first  half  of  the  cell,  where  it  strengthens 
the  shock  and  thus  helps  to  sustain  the  detonation.  Such  a mechanism  seems 
to  be  a crucial  one  in  order  that  the  detonation  can  be  self-sustaining. 

There  are  a number  of  requirements  given  here:  (1)  There  must  be  a collision 
of  transverse  waves.  (2)  The  hot  spot  does  not  arise  until  the  chemical 
reactions  become  exothermic.  (3)  The  resulting  pressure  pulse  must  reinforce 
the  lead  shock.  (4)  The  reinforcement  must  occur  in  the  first  half  of  the 
cell/73 

Barthe/83  has  shown  that  the  cell  width  or  spacing  correlates  rather 
well  with  pressure  at  various  dilutions  using  an  acaustic  theory  for  the 
origin  of  the  hot  spots.  The  theory  requires  that  caustics  associated  with 
two  rays  collide  at  the  hot  spots.  The  two  rays  start  out  in  different 
directions,  one  going  first  deeper  into  the  exothermic  zone  where  it  turns 
around,  then  travels  to  the  shock  where  it  is  reflected  and  subsequently  re- 
turns to  the  starting  plane;  the  other  first  starts  toward  the  shock,  re- 
flects, moves  past  the  starting  plane  and  deeper  in  the  exothermic  zone  to 
turn  around  and  return  to  the  starting  plane.  The  time  to  complete  this 
cycle  defines  a characteristic  time  t which  is  intimately  related  to  the 
characteristic  reaction  kinetic  times. 

Ill  INITIATION 

As  indicated  earlier,  detonations  are  usually  initiated  by  means  of 
igniters.  Clearly,  one  of  the  functions  of  an  igniter  is  to  generate  a shock 
wave  of  sufficient  intensity  so  that  the  reactions  become  fast  enough  to 
couple  rather  closely  to  the  shock.  If  one  can  ignore  the  effects  of  the 
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reactions  and  assume  the  igniter  energy  is  deposited  in  an  extremely  small 

fa} 

time,  then  one  can  use  point  source  blast  wave  theory v ' to  consider  the 
behavior  of  the  shock.  From  this  theory  R * t2/5  and  R “ t"3/5,  so  the 
shock  starts  out  at  a small  radius  at  an  extremely  high  velocity  but  sub- 
sequently is  always  decaying.  The  temperature- time  behavior  of  such  a wave 
is  shown  in  Fig.  1.  For  early  times  the  temperature  is  very  high  and  drops 
rather  rapidly  as  one  follows  the  particle.  With  reaction,  the  high  tem- 
peratures for  early  times  would  probably  require  that  the  reactions  would  al- 
ways be  endothermic  and  the  dropping  temperature  would  tend  to  slow  or  quench 
the  reactions.  When  the  shock  has  slowed  enough  so  that  the  reactions  be- 
come exothermic  at  the  end  of  the  induction  zone,  we  can  introduce  the  char- 
acteristic time,  T^nj.  For  the  pressure  pulse  associated  with  the  end  of 
the  induction  zone  to  propagate  to  the  lead  shock,  we  introduce  the  char- 
acteristic time,  T . During  the  time  x.  . + T the  lead  shock  must  not  have 
p ina  p 

decayed  too  much  if  a cyclic  process  is  to  be  started.  This  requirement 

means  then  that  there  must  be  a minimum  value  of  t for  the  blast  wave  of 

o 

Fig.  1,  which  in  turn  implies  a minimum  value  of  E,  the  total  energy  added 
to  create  the  blast  wave. 

The  studies  of  Bach,  Knystantas  and  Lee^10^  of  the  direct  initiation  of 
detonation  using  laser  pulsed  sparks,  ordinary  electrical  sparks,  exploding 
hot  wires  and  exploding  high  explosive  charges  as  igniters,  have  shown  that 
the  initiation  energy  required  for  direct  initiation  is  not  constant  from 
technique  to  technique.  Instead,  it  increases  both  as  the  duration  of  the 
time  for  energy  deposition  increases  and  as  the  kernel  volume  incxeases.  In 
fact,  they  found  that  the  important  parameter  for  energy  release  times  which 
span  the  20  nanosecond  to  60  microsecond  range  is  not  the  total  energy 


Fig.  1.  Temperature  along  pi 
Characteristics  and 
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deposited  but  the  rate  of  deposition  of  energy  per  unit  volume  of  the  ker- 
nel, i.e.,  the  power  density  (energy  per  unit  time  per  unit  volume)  of  the 
source.  Such  a result  seems  to  be  more  consistent  with  an  expanding  piston 
model  in  which  the  piston  drives  a shock  of  minimum  strength  and  the  power 
density  balances  the  rate  at  which  energy  is  expended  as  work  by  the  piston 
driving  the  shock.  In  this  viewpoint,  a maximum  shock  Mach  number  for  the 
lead  shock  is  reached  and  this  result  is  much  more  satisfying  than  starting 
off  with  an  infinite  Mach  number  as  one  must  do  when  one  uses  ideal  point 
source  blast  wave  theory.  This  experimental  result  further  suggests  that  the 
lead  shock  wave  cannot  be  modeled  as  a point  source  blast  wave  even  though 
the  energy  deposition  times  would  seem  to  indicate  that  it  should  be  possible 
to  model  the  shock  in  this  way. 

Bach,  Knystantas  and  lee  also  found  in  their  studies  that  if  the 
spherical  shock  wave  is  to  become  a detonation  wave,  then  transverse  struc- 
ture must  appear  on  a lead  shock  wave  that  previously  had  been  a smooth 
spherical  wave.  If  the  structure  did  not  appear,  the  wave  decayed  until  re- 
action ceased  and  transition  to  detonation  did  not  occur.  Such  a result  is 
quite  interesting  since  the  behavior  that  is  observed  when  initiation  is 
effected  is  similar  to  the  behavior  that  is  observed  during  the  propagation 
of  actual  self-sustaining  detonations.  It  thus  strengthens  the  notion  that 
the  transverse  structure  is  essential  for  the  propagation  of  self-sustaining 
detonations. 

The  appearance  of  structure  during  the  transition  to  detonation  again 
requires  consideration  of  Tg  if  an  acoustic  origin  for  the  hot  spot  structure 
is  correct.  The  value  of  r&  will  depend  upon  x^,  since  the  latter 
establishes  the  distances  over  which  the  sound  rays  must  travel.  Further,  for 


a constant  Mach  number  shock,  t is  greater  than  T , since  the  rays  for 

a p 

finding  the  value  of  the  former  are  inclined  to  the  flow  direction  as  well 
as  making  a complete  cycle  while  the  ray  for  finding  the  value  of  the  latter 
moves  directly  along  the  particle  path  toward  the  shock.  If  the  same  re- 
lationship holds  true  for  a decaying  wave,  then  the  parameter  (t^n(j  + Ta)/tQ 
is  a critical  parameter  for  the  establishment  of  structure  and  for  success- 
ful transition  to  detonation. 

IV  RESEARCH  PROGRAM 

There  are  four  distinctly  different  portions  of  the  research  program 
which  were  conducted  under  this  grant.  Two  of  these  are  experimental  and  two 
are  theoretical.  The  first  experimental  portion  was  planned  at  the  start  of 

I 

the  grant  and  consisted  of  a study  to  measure  the  structural  size  of  the  trans 
verse  wave  pattern  of  self-sustaining  detonations  in  mixtures  that  were  ap- 
propriate to  the  Air  Forces' use  in  FAE  devices.  This  is  reported  in  Section  V 
The  second  experimental  study  was  initiated  about  two  years  after  the  start  of 
the  grant  and  consisted  of  a reflected  shock  initiation  study  to  determine 
the  effect  of  promoters  and  inhibitors  on  the  delay  to  detonation  and  rate  of 
heat  release  in  propane-air  mixtures  and  also  to  determine  appropriate  delays 
to  detonation  and  heat  release  rate  for  propylene  oxide  and  oxygen-nitrogen 
mixtures  of  interest  to  Dr.  John  Lee  at  McGill  University.  The  results  of 
these  studies  are  reported  in  Section  VI. 

The  first  theoretical  effort  was  also  anticipated  at  the  start  of  the 
grant.  In  this  effort  an  initiation  pulse  (a  non- ideal  blast  wave)  was  to  be 
placed  in  a reactive  system  of  hydrogen,  oxygen  and  argon.  The  flow  fields 


were  then  to  be  computed  using  a full  set  of  reactive  equations.  On  these 
background  fields  an  acoustic  pulse  was  to  be  placed  and  subsequent  behavior 
was  then  to  be  followed  to  see  if  caustic  collisions  occurred  which  de- 
pended upon  the  position  and  timing  of  the  initial  pulse.  This  material  is 
described  in  Section  VII.  After  the  start  of  the  grant,  another  theoretical 
program  was  initiated.  In  this  program,  nonreactive  blast  wave  fields  pro- 
duced by  various  types  of  non-ideal  sources  were  used  to  monitor  the  initia- 
tion bohavior  using  hydrogen-oxygen  kinetics.  This  is  reported  in  Section 
VIII.  Additionally,  a number  of  studies  have  been  made  relative  to  the  nature 
of  non-ideal  blast  waves.  The  major  portion  of  these  studies  were  supported 
by  a grant  from  NASA  Lewis  Laboratories  and  the  Federal  Republic  of  Germany. 
However,  the  material  has  interest  to  the  Air  Force  needs  relative  to  ex- 
plosion behaviors  and  a summary  of  this  work  is  reported  in  Section  IX. 

V THE  EFFECT  OF  MOLECULAR  STRUCTURE  ON  TRANSVERSE  WAVE  SPACING  AND 
PROPAGATING  DETONATION 

As  an  initial  part  of  the  program  it  was  decided  to  study  the  effect 
of  structure  of  isomers  (same  molecular  formula,  different  structure)  on  the 
spacing  of  transverse  waves  for  the  isomers  of  ethylene  oxide  and  propylene 
oxide.  Tables  1,  2 and  3 list  the  thermochemical  properties  and  structure  of 
these  molecules.  Ethylene  oxide  has  only  one  isomer,  acetaldehyde,  however 
propylene  oxide  has  five  isomers  which  are  listed  in  Tables  1,  2 and  3.  Also 
in  Table  1 the  abbreviations  that  we  used  for  these  different  species  are 
given. 

The  results  of  this  study  were  not  very  definitive.  Figures  2 and  3 
ropresent  plots  of  the  cell  spacing  measured  for  these  various  species  as  a 
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TABLE  3 


MOLECULAR  STRUCTURES  AND  PERTINENT  DATA 
OF  C3H60  ISOMERS:  mw  58.078 


PROPYLENE  OXIDE: 

-6-6-6 

b' 

TRIMETHYLENE  OXIDE: 

V-°^c' 

* \ 

ACETONE: 

-c-c-c- 

6 

PROPIONALDEHYDE: 

I It 

-c-c-c=o 

I I 


ALLYL  ALCOHOL: 

-C=C-C-OH 


V.P  ■ 555  torr 

B.P  a 35° C 

AHf  (g)  *'22.17  kcal/mole 

V.P  « 325  torr 

B.P  = 47  • C 

AHf  (Q)  * - 19.25  kcal/mole 

V.P  =229.4  torr 

B.P  = 56.5  °C 

AHf  (g)  = -52.00  kcal/mole 

V.P  =285  torr 

B.P  =48.8  °C 

AHf  (g)  =-45.90  kcal/mole 

V.P  =28.1  torr 

B.P  =97  °C 

AHf  (g)  =-31.55  kcal/mole 


METHYL  VINYL  ETHER: 

-6-0-6= C 

i 


V.P  =2208  torr 
B.P  = 8.0  # C 
AHf(g)  = — 


ALL  VAPOR  PRESSURES  (V.P)  ARE  AT  25  *C 
ALL  AHf  ARE  AT  298 °K 


2 


length  versus 


function  of  both  the  final  shock  temperature,  i.e.,  the  Chapman-Jouguet  tem- 
perature of  the  detonation, and  as  a function  of  the  initial  temperature  be- 
hind the  shock,  i.e.,  tho  nonreactive  shock  temperature.  Here,  the  tem- 
perature was  varied  by  varying  the  stoichiometry  and  dilution  in  the  experi- 
ments. The  most  significant  and  striking  differences  are  the  fact  that  the 
slope  of  the  spacing-shock  temperature  curves  is  extremely  high  for  both 
ethylene  oxide  and  propylene  oxide  and  relatively  flat  for  all  of  the  isomers. 
The  conclusion  that  can  be  drawn  from  this  study  is  that  both  propylene  oxide 
and  ethylene  oxide  are  the  most  reactive  isomers  in  their  series.  Details  of 
this  work  are  reported  in  a Master's  thesis  written  by  Mr.  Delaney. 
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VI  REFLECTED  SHOCK  INITIATION  STUDIES 

Two  distinct  studies  were  performed  in  this  part  of  the  program. 

These  were  initiated  by  Mr.  Joe  Foster  at  Eglin  Air  Force  Base  approximately 
two  years  after  the  program  started.  The  first  of  these  was  a study  of  the 
effect  of  inhibitors  or  promoters  on  the  ignition  delay  in  a propane-air  mix- 
ture behind  a reflected  shock  in  a conventional  shock  tube.  The  second  of 
these  was  a similar  study  of  initiation  behind  the  reflected  shock  for  a 
mixture  of  propylene  oxide,  oxygen  and  nitrogen. 

The  experimental  work  for  both  of  these  studies  has  been  completed. 
Eventually  these  two  studies  will  yield  two  Master's  theses  describing  the 
results  of  the  studies.  These  theses  are  not  completed  at  this  time,  but 
should  be  completed  by  the  Fall  of  1977.  Copies  of  these  theses  will  be 
supplied  to  the  Technical  Monitor  at  the  time  when  they  are  completed. 
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VII  ACOUSTIC  WAVES  ON  A REACTIVE  FLOW 


Earlier,  we  noted  that  if  a decaying  spherical  shock  wave  is  to  be- 
come a detonation  wave,  transverse  structure  must  appear  on  the  initially 

(O') 

smooth  shock.  Also,  Barthel^  J had  proposed  an  acoustic  mechanism  whereby 
colliding  caustics  formed  hot  spots  tc  initiate  the  transverse  structure  in 
steady  planar  detonations.  Here  we  examine  whether  a similar  mechanism 
might  be  possible  for  the  spherical  case. 

7.1  Ray  Equations 

This  approach  requires  that  we  follow  various  rays.  Hence,  following 
f 121 

Milne,  J we  write  the  ray  equations  as 

, 5^  ’ ani  + ui  i - 1,  2,  3 (1) 


n. 

i 


Rnl  3 a 
dt  3x. 


+ n 


3jV 

k 3 x. 


pJLl 

dt 


3a 

3x. 


+ n. 


3uk 


(2) 


where  a is  the  local  speed  of  sound,  n^  is  the  component  of  the  unit  normal 
to  the  wavefront  and  ui  is  the  velocity  component.  Upon  transforming  to 
the  spherical  case,  these  become 


u + ail 


(3) 


if 


am/r 


dS, 

dt 


♦ l 


3 u 

57" 


(a+£u)/r 


(4) 

(5) 


where  l and  m are  the  components  in  the  r and  ^ directions  of  the  normal  to 
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the  wave  front  and  u is  the  particle  velocity  in  the  radial  direction. 

Now  consider  the  behavior  of  (a+Jtu)/mr  as  a function  of  time  along  the 
ray  path.  Differentiating  the  expression  as  we  follow  the  ray  path,  we  get 


i_ " atfa,*u)  - ft (mr) 

dt(  mr  j m4r4 


„/3a  A 

3 a 

dr  , 

o 3u 

. 3u  dr 

x dA  ) 

x o..)  L dm  . 

dr) 

- 

H 

3t 

*St 

fr  dt] 

+ u dtf 

- a + £u  r + 

X — 1 ■ — 

m dt 
— ...  w- 

"nrr7 


(6) 


0 J _ 

Using  * ’ in  3t  an<*  su**stituting  Eqs.  3,  4 and  5 into  Eq.  6,  we  get 


d a+&u|  m 1 
dt  mr  J " mr 

Alternatively,  the  left-hand  side  of  Eq.  7 can  be  written 


da 

5t 


+ A 


du 

Si' 


j 

t * 


(7) 


a+Au) 

_1  d fa»Au 

a+Au  dr 

mr  J 

r dt  ( m t 

mr  dt 

Combining  Eqs.  7 and  8,  then  multiplying  by  r,  we  get 

d fa+fcu)  _ (a+Au)  dr  „ 1 [da  , |u) 
3t  l m J mr  cEt  m |Tt  * 5Tj 


(8) 


(9) 


Hence,  if  r is  large  and  the  right-hand  side  of  Eq.  9 is  zero,  we  obtain  the 
planar  equation  for  steady  flow 


(10) 


which  has  the  solution 


Si**.  . c 

m 


(ID 


where  c is  a constant  for  each  ray. 
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If  the  spherical  a,  u fields  are  similar  to  those  found  in  a steady 
planar  detonation,  then  we  can  expect  colliding  caustics.  Therefore,  for 
large  r,  we  expect  structure  in  the  spherical  case. 

The  behavior  for  small  r is  not  so  apparent,  since  we  have  to  solve 
Eq.  7 and  it  contains  unsteady  terms.  If  these  unsteady  terms  are  so  small 
that  they  are  negligible  compared  to  r,  we  obtain 


a + Jtu 

TT"  ci 


02) 


where  c^  is  a constant  which  is  an  angular  rate.  For  perhaps  some  par- 
ticular a,  u fields,  we  believe  this  equation  could  yield  colliding  caustics. 
However,  whether  the  unsteady  fields  created  by  rapid  energy  addition  in  a 
reactive  gas  would  create  such  a field  can  only  be  determined  by  numerical 
integration  because  of  the  complexity  of  the  equations. 

Our  previous  success  using  an  implicit  integration  technique  for  the 
steady  planar  detonation  case  led  us  to  believe  that  we  might  be  able  to 
adapt  our  previous  program  to  the  unsteady  spherical  case  with  energy  ad- 
dition in  the  center  of  the  reactive  mixture. 


7.2  Equations  for  Reactive  Flow 

he  assumed  no  diffusion  or  heat  condition  occurred  in  the  mixture. 
Then  the  species  equations  for  a reactive  flow  are  given  by 

D[AJ  DIAj), 

“DiT-  * -Dri*  * lAiJlur + j ’ 1-1 NSP  C13) 

OIAJ, 

where  JA^.]  is  the  molar  concentration  of  species  i,  — L is  the  rate  of 
change  of  concentration  due  to  chemical  reaction  as  we  follow  a fixed  mass 


ff 
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of  mixture,  j » 0,1,2  for  planar,  cylindrical  and  spherical  cases,  u and 


ur  are  the  particle  velocity  and  its  spatial  derivative  and  r is  the  position 


D d 3 

coordinate.  Also  j5t  " St*  + u STr  *s  t^'e  ®uler  derivative  as  we  follow  the 


mass. 


For  our  system,  the  perfect  gas  law  takes  the  form 
NSP 


P - l [A.]  RT , 
i-1  1 


(14) 


where  p is  the  pressure,  R is  the  universal  gas  constant  and  T is  the 
absolute  temperature. 

The  molar  and  mass  densities  are  given  by 


NSP 

"n  " J. 


(15) 


i=l 


p _ N?P  MW.JA.J  , 


l 

i-1 


(16) 


where  MW^  is  the  molecular  weight  of  species  i. 


From  Eqs.  13  and  15,  p^  obeys  the  equation 


DP  NSP  D[A  ] 

I "StH  ' °ntur  * > F> 

i-1  ut  R n r r 


Dt 


(17) 


while  from  Eqs.  13  and  16,  p obeys  the  equation  for  conservation  of  overall 


D£ 


Us 


Dt  “ “ p(ur  + * * 


(18) 
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The  first  law  is  written  as 


Dh 

Dt 


I DR 

p Dt 


(19) 


NSP 


whore  h ■ l (A.Jh./p  is  the  enthalpy  per  unit  mass  of  mixture,  h,  is  the 
i-1  1 1 1 

enthalpy  per  mole  of  species  i,  and  q is  the  rate  of  energy  per  unit  mass 


added  from  outside  the  system.  We  further  assume  that  h,  is  a function  of 


3^ 


i 

T only,  with  hj  « hjfT)  and  ipp  ■ Cp  (T). 


From  Eq.  19  and  the  relationship  for  h,  we  obtain 


DT 

Dt 


NSP 
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i=l 


hi 

T-XI 
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NSP 
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The  momentum  equation  is  given  by 


Du 

Dt 


P 


(21) 


The  particle  path  is  given  by 


Dr 

Dt 


(22) 


The  right-hand  side  of  the  equations  for  [A^] , T,  u and  r contain  the 
variables  [Aj],  T,  u,  r,  uy  and  pr  if  q is  specified  as  a function  of  these 


variables  as  well  as  of  time.  For  explicit  programs,  ur  and  pf  need  only  be 


calculated  at  time  t ■ t to  find  the  variables  at  t » t .,  for  we  have 

n n +1 


NSP  + 3 equations  to  calculate  NSP  + 3 dependent  variables.  However,  the  set 


i 

* 


of  equations  is  mathematically  a set  of  "stiff"  differential  equations  which 
often  exhibit  numerical  stability  problems  during  their  integration  unless 
very  small  step  sizes  are  taken.  In  an  effort  to  avoid  this  difficulty,  we 
turned  to  an  implicit  method. 


7.3  Implicit  Methods 

The  concept  of  the  implicit  method is  basically  as  follows.  We 
have  a set  of  differential  equations 

dyi 

- fiC/j.t)  i,j  - 1,  , N . (23) 


We  know  the  values  y.  at  time  t » t . We  wish  their  values  at  time 

J n 

t . Vl.  bet  us  define  " Vl  ‘ ln  ,nd  "i.n.l  * ’'i.ntl  ‘ yi,n'  Then- 
we  assume  an  implicit  expansion,  such  as 


Ki,n+1 


1 

[dyi.»*l  , dyi.n' 

h . d V1 

+ 1 

'dyi,n.l 

dy. 

2 

l dt  dt  J 

n+1  dt 

2i 

dt 

dt 

(24) 


24  is  implicit  because 
, so  are  not  known  yet. 


dyl.n.l 

dt 

But 


depends  upon  the  values  of  y 
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at  time 


dyi.n*l 

dt 


dy 


dt 
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i.n  . 
dtz  hn+l 
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jjnjl  dyi,n  „ „ . * V a W7J.n  . 

dt  ' gi,n  n+1  dt  n+1 
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dy. 


dt 


j-1 


(26) 


If  we  use  Eq.  24,  Eq.  26  becomes 


dyi,n+l  dyi,n  . . 8 . L hn»lp 
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N 


where 


a,  , « l 0.  .0.  . . Putting  Eq.  27  back  into  Eq.  24,  we  get 
K»J 


i 

i .}< 


I } 


s .a 


i,n+l  2 


h«*i  N 

n+1  j 


j=l 


6iJ 
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2 i,j 
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J >n 


, h ,k 

^i,n  n+1  2 


h , N 
n+1 


gi,n  2 .^0i,jgj,n 


(28) 


which  is  a coupled  set  of  N linear  algebraic  equations  in  Kj  n+^.  Upon 
solving  for  n+1,  then  yA  n+1  is  easily  obtained  from 


yi,n+l  “ yi,n  + Ki,n+1  * 


(29) 
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7.4  Closure  Considerations 


Du 


Dp, 


Prom  the  preceding  theory,  it  appears  that  equations  for  and 


are  needed.  We  attempt  now  to  derive  them  by  using  the  relationship 
DYr 


.3  D£ 

Dt  1r  Dt  Ur^r 


(30) 


Applying  this  equation  with  $ « u,  we  get 


Du 


Dt 


3 

r~ 

4r 


u‘ 


P--  P„P- 
SL  + .r...r-  - u2 
p pz  r 


(31) 


Thus  we  now  have  two  additional  variables  p_  and  p . Clearly,  to  follow 

rr  r 


such  a scheme  will  not  lead  to  closure. 

Rather,  closure  is  to  be  accomplished  by  a different  scheme.  We 

I 

integrate  along  a path  from  the  point  i-l,n  to  the  point  i,n+l  or  from 
point  i+l,n  to  the  point  i,n+l  by  using 


dKr 


It 


* u + K+a 


(32) 


and 


6Kr 


6t 


* u - K'a 


(33) 


which  is  much  like  integrating  along  characteristics  except  it  is  more 


convenient  to  move  from  computed  points  on  the  t **  tR  time  line  to  avoid 


interpolation  between  points  on  this  time  line  as  required  when  integrating 
along  characteristics.  Further,  because  we  now  have  two  new  variables  K+ 
and  K~,  we  also  integrate  the  equations 
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and  a in  Eqs.  32-38  is  the  frozen  speed  of  sound  given  by 

c 


2 i=l 
a = 


NSP  p 

2 IAi]  nr 
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JjA  iJ 


i*l 


Pi 
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p 

-&RT 

P 


(39) 


Now  it  appears  that  the  system  is  overspecified,  but  the  extra  equations  are 
used  for  checking  purposes.  However,  we  also  will  note  here  that  care  must  be 
exercised  in  the  integration  of  Eqs.  32-38  but  the  discussion  of  this  problem 
will  be  deferred  to  a later  section.  We  see  though,  that  upon  integration,  one 
can  find  expressions  for  K",  K+,  u (i,n*l)  and  p_(i»n+l)  in  terms  of  {A.],  u, 
r and  T at  the  new  time,  so  sufficient  equations  are  at  hand. 

7.5  Mixture  and  Reactions 

For  all  calculations,  the  system  was  assumed  to  be  70%  argon  and  30% 
of  a stoichiometric  mixture  of  hydrogen  and  oxygen.  For  this  system,  the 
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raactions  and  rate  constants  given  in  Table  4 were  used.^*^ 

7.6  Initial  Conditions  and  Boundary  Conditions 

The  initial  conditions  were  that  the  mixture  was  at  rest,  with  uniform 
composition,  pressure  and  temperature  throughout.  The  initial  pressure  was 
200  Torr  and  the  initial  temperature  was  300*K. 

The  boundary  conditions  were  specified  at  r ■ 0 and  at  the  moving  wave 
forming  the  outer  boundary. 

At  r * 0,  symmetry  requires  u ■ pr  * fA^]r  « Tf  ■ 0. 

The  outer  boundary  may  be  either  an  acoustic  wave  or  a shock  wave. 

If  the  former,  its  path  is  given  by 


rf  * rc(0)  + aot  , (40) 

where  rf  is  the  position  of  the  front  at  time  t and  rc(0)  is  the  origin  of 
the  wave.  Because  of  the  nonlinear  nature  of  the  flow,  the  acoustic  wave 
will  gradually  get  steeper  and  ultimately  become  a shock.  To  determine  when 
this  transition  occurs,  we  assume  no  reaction  and  note  that  before  transi- 
tion, p£  ■ pQ  and  u£  * 0 for  all  points  along  the  front.  Whence 


S+Ui 

ir\ 


f 


p^"  + aour£  " 0 ' 


(41) 


so  p *»  o a u (42) 

Frf  Mo  o rf  K J 

on  the  front.  Since  Eq.  42  holds  for  every  point  on  the  front,  then 


d*pri  6+u 

»' ■■■  * pa_  r 

Ot  |£  OO  — 


(43) 
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TABLE  4.  The  reaction  equations  and  rate  constants  used 
in  tho  Rayleigh  program.  The  rate  equations 
are  of  the  form  ■ Z^Tn  exp  (E^/RT) . 


Reactants 

Products 

2ia 

n 

Eil 

02  + M 

0 ♦ 0 ♦ M 

3.56  D 15 

0 

1.18  D 05 

0 + 0 + M 

02  ♦ M 

2.00  D 14 

0 

0.0 

R,  ♦ M 

H + H + M 

3.10  0 IS 

0 

1.10  O 05 

H + H + M 

h2  + M 

6.50  D 14 

0 

0.0 

°2  + H 

0 + OH 

2.24  D 14 

0 

1.68  D 04 

OH  + 0 

°2  + « 

1.30  D 13 

0 

0.0 

H02  + M 

02  + H + M 

2.40  D 15 

0 

4.59  D 04 

02  ♦ H + M 

H02  + M 

1.59  0 15 

0 

-1.00  D 03 

h2  ♦ ° 

OH  + H 

1.74  D 13 

0 

9.45  D 03 

H + OH 

h2  ♦ 0 

7.33  D 12 

0 

7.30  D 03 

H20  + H 

H2  + OH 

8.41  1)  13 

0 

2.01  D 04 

H2  + OH 

HjO  + H 

2.19  D 13 

0 

5.15  D 03 

h2o  +'  0 

OH  + OH 

5.75  D 13 

0 

1.80  D 04 

OH  + OH 

h2°  ♦ o 

5.75  D 12 

0 

7.80  D 02 

H20  + M 

H + OH  + M 

1.00  D 17 

0 

1.17  D 05 

OH  + H + H20 

h2o  ♦ h2o 

1.17  D 17 

0 

0.0 

OH  + H + H2 

«2°  + H2 

2.92  D 16 

0 

0.0 

OH  + H ♦ 02 

h20  ♦ 02 

2.92  D 16 

0 

0.0 

OH  + H + AR 

H20  + AR 

7.02  D 15 

0 

0.0 

H2°2  + H 

H2  + ,,02 

2.34  D 13 

0 

9.20  D 03 

h2  ♦ ho2 

H2°2  + H 

9.60  D 12 

0 

2.40  D 04 

H202  + OH 

h20  + H02 

1.00  D 13 

0 

1.80  D 03 

H20  + H02 

H202  + OH 

2.80  D 13 

0 

3.27  D 04 

H202  ♦ M 

OH  + OH  + M 

1.20  D 17 

0 

4.55  D 04 

OH  + OH  + M 

H202  + M 

8.40  D 14 

0 

-5.30  D 03 

H2°2  + H 

H20  + OH 

3.18  D 14 

0 

9.00  D 03 

H20  + OH 

H2°2  + H 

5.60  D 13 

0 

7.79  D 04 

a(cm*/mole/sec  or  cm6/mole2/sec) . 
k(cal/mole) . 
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along  tha  front.  Consequently, 


Trl,  ■ 4-  (§?}/  Vr|/  “P: 


f 


rr!f  " " 


M?).' 


YPr 


U +j 

rr 


u 


r u 
r “ r 


- YPr  Vlq  ' ur(Wr)  4 Vrr,  “ B?  * “J  4 *V 

f f r L 


Du  2 
- u * + 


■]. 


p a 
Ho  o 


Prr  Pr  pr 
**£  r rf  a 


1 u;  + a u 


o rr. 


(44) 


We  have  also  made  the  assumption  that  q * 0 at  the  front,  so  the  flow  is 


Pr  YPr 

isentropic  in  that  region  and  therefore  — » -—  . Consequently 

K P„u.  P P 

rf  0 rf 

p * rr~  3 — - — leads  to  the  following  relationship  between  p and  u 
rf  ao  ao  rrf  rrf 


from  Eq.  44. 


P*  p a 
rr^  Mo  o 


u 


rf  , rf 
u__  + - — + *r  - — 


rrf  ao  2 rf 


(45) 


6+u 


Putting  this  result  into  the  equation  for  — , we  have 
6+u 


~xrfm 


u 


2 + i 


a. 


r-  2 o r 


(46) 


which  can  be  written  as 


6+ 


) 


. i!a 

2rf 


i 


U- 

lrfj 


m 

2 


(47) 
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Eq.  47  has  the  integrating  factor  exp 
so  the  solution  of  Eq.  46  is 

♦ lii  f _ dt 

V'O 


if^dt 
2 J 


- [rc(o)  + aQt] 


-i/2 


1 


ur[»c(o)*a„t]j/*  ’ Cl  * » J [Wli/2  ' 


(48) 


For  the  spherical  case  with  j ■ 2,  the  last  integral  becomes 
An[rc(o)+aQt].  If  at  t * tj.  uy^  - u^. 


u 


rcOO+V 


2 

Y+l 


2 

JCil 


1 


V(rc(o)Vl)  +4'o  £n 


rcCo)+V 


rcf°)+Vl 


(49) 


If  r « 0,  this  becomes 
c 

2 

y+i 


rf 


r 2 

Vi 


(50) 


But  the  shock  wave  forms  as  u ■+■«>,  or  from  Eq.  50, 

rf 


exp 


JL 

JCi. 


V» 


(51) 


and  since  u < 0,  a shock  forms  in  a finite  time  larger  than  t,. 
rl 

Alternatively,  it  may  happen  that  before  this  value  of  t is  reached, 
an  interior  shock  overtakes  the  leading  acoustic  front.  In  that  case,  there- 
after the  shock  equations  must  be  satisfied  at  the  outer  boundary,  with  con- 
stant pressure,  density  and  temperature  and  zero  velocity  ahead  of  the  shock. 
Since  these  conditions  pertain  to  a special  case  of  a shock  wave  moving  into 


r Vifiitf 
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a reacting  fluid,  the  theory  for  this  general  case  is  considered  next. 

7.7  Shock  Wave  in  a Reacting  Fluid 

In  a system  moving  with  the  shock,  the  conservation  equations  may  be 

written 

paws  " pbIws  4 <Vub)]  (52) 

Pa  + paws  ' pb  + pbIws  * <Vub)]2  (53) 

ha  + </2  * \ * Iws  * (ua'ub)]2/2  , (S4) 

where  the  subscripts  a and  b refer  to  conditions  ahead  of  and  behind  the 
shock,  wg  is  the  velocity  of  the  shock  relative  to  the  fluid  ahead  and  the 
plus-minus  signs  refer  to  outwardly  and  inwardly  moving  shocks.  We  also 
have 


h * 

NSP 

l 

i«l 

[A.l^/p 

(55) 

P • 

NSP 

l 

i«l 

[A.]MW. 

(56) 

Pn 

NSP 

l 

i-1 

[Ai] 

(57) 

P - 

p RT 
n 

(58) 

»f  ■4*r 


(59) 
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NSP 

l ‘Vv 


c « 
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i»l 


cv 
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cp  - cv  - R 
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We  further  make  the  assumption  that  the  flow  is  frozen  while 
crossing  the  shock,  so 


Hence 


[Vb  tA,], 


“b  0. 


and 


\ \ 


pb  pa 


(60) 


(61) 


(62) 


(63) 


(64) 


(65) 


(66) 


which  means  that  the  molecular  weight  does  not  change  on  crossing  the  shock. 
There  are  now  sufficient  equations  for  solving  for  the  flow  conditions  behind 


the  shock  in  terms  of  w„  and  flow  conditions  ahead. 

s 


In  order  to  find  the  spatial  derivatives  uy  and  pf  behind  the  shock, 


we  differentiate  Eqs.  52,  53,  54,  55,  58,  65  and  66  while  following  the  shock 
to  get 
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fiV  6Sw  fiSp.  fiSw  6Su  6su, ' 

tt “s  * pa  Tt — rr  tv("a-"b»  *»bir‘  -5 — sr. 


6sp  fisp  fisw  fisp  6sp. 

•7HT  + “ST  WS  + 2paws  -1ST  = TT  + ~5t“  Iws±cWJi 


6sw_  ffisu  6su. 


' C Q Kl 

2 Pb[ws±Cua-ub)]  -j^-4  g^i 


fiSw  Ssh. 

S D 


F+Wsirsx+  [Vcvub)]  ir 4 rre 5t~j 


6Sws  . [6\  6\ 


f\  £tAi]cPi  «\  H - \ 

t p 6t  pb  it  ^ pb  St 


u ^b 

St~ 


o T.  n. 

R <\  IT  * Tb  Tt* 


1 6S[Ai^b  ^b  5 pb  1 ^Va  tAA]a  6SPS 


fit  Pa  fit 


■pf  St 


1 S pnb  pnb  {5pb  . 1 "a  "a  jSpa 

Pb  fit  p£  fit  pa  ot  p*  fit 


where  the  derivatives  ahead  of  the  shock  have  the  form 


6\  Dip  , Dip  U 
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derivativos  behind  the  shock  have  the  form 

Ss<|>.  Dik  a ^ . Dip,  dip, 

-r‘r‘('s-ub)')T!ri[,!!  (vV]  rr 

6 s w 

and  the  derivative  — g~  is  given  by 


(75) 


6sw 

~sr  “ ♦ 


i'.-^  ■ 


(76) 


with  r being  the  shock  position  coordinate. 

s » s 

6 pb 

The  set  of  (NSP+6)  equations,  67-73,  is  now  solved  for  — and 
6su.  6sw 

■■•g--"  in  terms  of  and  derivatives  ahead  of  the  shock.  When 

6sp  6su 

and  — are  written  in  the  form  of  Eq,  75,  they  contain  only  the  two 

t 

unknowns  u and  p , so  simultaneous  solution  of  these  two  equations  will 
rb  rb 

yield  values  for  u and  p . 

rb  rb 


7.8  Specification  of  q 

The  rate  of  energy  input  was  assumed  to  have  the  form 

q - qQ  X (r,rcCt))<f>(t)  , (77) 

where  qQ  is  a constant  and  rc(t)  is  the  outer  edge  of  the  core  in  which 
energy  is  being  added.  For  all  work  undertaken  in  these  studies,  4>  had  the 
form 

(t)  « exp(-t/tj)  - exp(t/t2)  , (78) 

where  Tj  and  t2  are  decay  parameters  or  shaping  constants.  In  addition, 

r (t)  had  the  form 
c 
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rc(t) 


V 


a 


(79) 


with  f being  a constant  and  the  exponent  ct  was  within  the  range  0 < a < 1, 

Initially,  X(r,r  ) was  chosen  to  be  the  function 

C 

X « r*  - rs  . 0 < r < r (80) 

c c 

* 0 r < r < « i 

c 

This  choice  for  X turned  out  to  be  an  unfortunate  one,  for  at  r » r there 

c 

is  a jump  in  qf.  This  discontinuity  led  to  difficulty  in  the  computing  pro- 
gram. The  specific  difficulty  and  the  theory  will  be  given  in  later  sections.  i 

i 


7.9  Additional  Computational  Matters 

Since  an  implicit  method  for  integration  of  a set  of  differential 
equations  is  usually  slower  and  more  costly  than  an  explicit  method,  to  keep 
costs  down  the  maximum  number  of  particles  to  be  followed  was  limited  to  60. 
When  the  number  of  particles  followed  approached  60,  the  flow  field  was  re- 
zoned to  a smaller  number  before  calculation  at  the  next  time  line  started. 

Four  types  of  points  occurred  in  the  calculation,  namely,  (1)  origin, 
(2)  interior,  (3)  shock  and  (4)  outer  boundary.  The  treatment  at  the  origin 
had  to  be  special  because  the  boundary  conditions  require  that  the  particle 
motion  and  momentum  equations  be  dropped  from  the  set  to  be  integrated. 

Shock  points  and  outer  boundary  points  also  had  to  be  treated  as  special 
cases. 


With  the  above  scheme  and  the  assumed  X function,  difficulty  in  start- 
ing the  computation  was  experienced.  Considerable  effort  was  devoted  to 
trying  to  use  similarity  variables  to  overcome  the  starting  problem,  but  this 
effort  was  unsuccessful.  Finally,  because  of  doubts  about  whether  the 
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implicit  method  could  be  used  for  starting  the  computation*  we  turned  to  the 
CLOUD  program  of  A.  K.  Oppenheim  as  modified  by  A.  A.  Adamczyk. After 
a few  time  lines  were  calculated  using  this  program,  negative  velocities 
began  to  appear  near  the  origin  and  unstable  oscillations  appeared  in  the 
flow.  These  oscillations  usually  became  so  large  that  the  program  would  fail. 
The  origin  of  this  behavior  was  traced  to  the  particular  X function  assumed. 
When  a new  X function  was  chosen,  given  by 

X ■ r * - 3rcra  + 2rs  , (81) 

these  troubles  disappeared.  Tt  should  be  noted  in  this  case  that  qy  is 

zero  at  r « r . 

c 

A second  difficulty  encountered  in  using  the  CLOUD  program  was  that 
u^  and  pr  derived  from  curve-fitting  the  values  of  u and  p as  a function  of 
radius  did  not  vary  smoothly  with  radius.  The  difficulty  was  traced  to 
roundoff  error  and  overcome  by  recasting  the  calculations  for  the  dimension- 
less pressure,  specific  energy  and  specific  volume,  p,  e and  v respectively, 
in  terms  of  the  reduced  variables  p-1,  e-e0  and  7-1,  where  eQ  was  the  initial 
dimensionless  specific  energy,  . 

One  of  these  time  lines  was  chosen  to  be  the  starting  line  for  the 
implicit  program.  The  particular  time  line  chosen  was  one  in  which  the  tem- 
perature in  the  core  had  been  raised  less  than  0.5®K,  so  no  significant  re- 
action could  have  yet  taken  place. 

7.10  Effects  of  Discontinuities  in  q and  q„ 

^ r 

This  section  summarizes  results  submitted  elsewhere.  ; Basically, 
we  sought  to  explore  what  effect  discontinuities  in  q or  qf  occurring  along 
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rc  • u+ka  would  have  on  the  solutions  of  the  set  of  equations 


ft  - vtur  * $ , 

(82) 

§t  ■ -VPr 

(83) 

§f  ■ -CY-l)eCur  + jp)  + <1  » 

(84) 

§£  - -YP(ur  ♦ j~)  + CY»l)q/v  , 

(85) 

pv  - (Y-l)e 

(86) 

near  the  discontinuity. 

The  representation  of  these  discontinuities  in  q or  ^ are  shown  in 
Fig.  '4.  The  particular  effects  on  p,  u,  e and  c are  dependent  upon  the 
value  of  k.  If  k * 0,  then  we  have  a behavior  like  a contact  discontinuity. 
If  k t 0,  then  the  effect  is  wavelike. 

A jump  in  a variable  will  be  indicated  by  the  usual  convention  of 
brackets,  such  as  ' ( 

[u]  - u2  - ux 

v 

where  the  subscripts  2 and  1 refer  to  points  just  inside  and  outside  the 

r » r surface, 
c 

The  derivative  following  a variable,  \p,  along  the  interface  is  given 


by 


If  ■ IS  * If 


(87) 


Also,  we  repeat  the  equation 


9 

W 


at  . 

Dt 


u tb 
ryr 


(88) 


Fig.  4.  The  assumed  behavior  of  the  energy  addition  rate  q vs  radial 
distance  r.  For  (a)  q and  q^  are  both  discontinuous  at  the 
kernel  edge  while  for  Cb)  q is  continuous  but  qy  is  discontinuous 
at  the  edge. 
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a)  Contact  Surfaco 

In  this  case,  k » 0,  so  r£  ■ u and  pressure  and  velocity  must  be 
continuous  across  the  interface.  For  a jump  in  q,  we  form 


^^•“YP2Iur3  + CY-Uq2v2  • 0 , 

[ur3  ■ (Y-l)q2/YP2r2  ■ q2/Y«2 


Using  this  result,  we  can  find  that  the  equation  for  {ej  is 


£|2l*  -CY‘1)  I®]  Cur^  + jux/rc)  ♦ q2/Y  , 


which  formally  has  the  solution 


[ej  * expi-ACtjHcj  + ^ 1 q2exp{A(T)}dT}, 


where 


ACt)  - (Y-l) 


j Cur  + U1 
jo  rl  1 


/rc)dt 


With  quiescent  initial  flow,  c^  is  zero  and  we  find  that  lej  is  positive, 
which  intuitively  is  the  expected  behavior. 


Prom  Eq.  86  we  get 


so  [v]  is  also  positive. 


[v]  - Vj/e^e] 


Turning  now  to  Eq.  83,  we  write 


"(v2Pr2  - v,^)  - 0 . 


Hence,  there  is  no  jump  in  vpf  but  there  is  a jump  in  py  given  by 


tPr1  ' (^-l)  % * - JrjM* 
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which  is  positive,  since p_  is  negative  and  [ej  is  positive. 

rl 

The  results  are  depicted  in  Pig.  5. 

For  q * 0 but  a jump  in  at  rc,  we  have  no  jumps  in  ur,  e,  v and  p^. 
We  now  examine  the  next  higher  derivatives  of  these  variables.  We  find  that 

fiIPJ 

-37—  " ‘YPl^rr^  * CY-l)qr2/v2  * 0 (97) 

so 

lUyyJ  ■ (Y-Dq,  /YPiVj  - q*  h ^ , (98) 

2 2 

which  is  negative,  since  q_  is  negative. 

2 

The  equation  for  Je^J  is  obtained  as 

“IT^ " -CY-l)Ier]{u  + J^/rJ  + ^ h , (99) 

1 1 2 

which  is  of  the  same  form  as  Eq.  91.  It  follows  then,  that  if  [erJ  ■ 0 at 

t **  0,  then  the  sign  of  Je  J depends  upon  the  sign  of  q*  • Since  the 

r t2 

latter  is  negative,  then  X«rJ  is  negative. 

For  Iv  ) we  get 

Ivr]  - vif®r^ei  ' C100) 

which  is  also  negative. 

Finally,  forming 

filuj 

TT"  ‘ -,ri[,r]  .0  , £101) 

we  get 

XPrrJ  - “Pr^vrJ/vl  “ “Pri^*r]/*1  * C102) 

which  is  negative. 

These  results  appear  in  Fig.  6. 


0 


0 


(e) 

Fig.  5. 


Property  variations  as  a function  of  radiu*  when  both  q and  qy 
are  discontinuous  at  the  kernel  edge  and  that  edge  moves  along 
a particle  path,  (a)  uy  vs  r,  (b)  py  vs  r,  (c)  e vs  r and 
(d)  v vs  r. 
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Fig.  6*  Property  variations  as  a function  of  radius  when  q is  continue 

but  is  discontinuous  at  the  kernel  edge  and  that  edge  moves 

along  a particle  path,  (a)  uf  vs  r,  (b)  pr  vs  r,  £c)  «r  v*  r i 

(d)  v.  vs  r. 
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b)  Wave  Surface 

For  k 0,  so  long  as  there  is  no  shock  wave*  u,  p,  v and  e are 

tinuous  across  rc,  with  no  jumps  and  their  derivatives  are  aero. 

For  a q jump  at  r . from  Eqs.  83,  85  and  87,  we  find 
c 

-Vjtpr]  - (rc  - upiur3  - 0 


- -YP2Iur3  + CY-l)q2/vx  - (rc-u1)IprJ  - 0 . 


Simultaneous  solution  then  yields 


lur]  « CY-lDqj/CaJ-C^-Uj)2)  - (y-l)q2/a*(l-ka)  , 

which  is  positive  for  k < 1 and  negative  for  k > 1,  and 

[Pr]  - (Y-llkc^/VjajCl-k2)  , 

which  has  a similar  behavior. 

From  Eqs.  84  and  87,  we  find 

-CY-DejCUy]  + q2  + (rc-Uj)fer]  » 0 , 

so 

I«r3  - -q^Y"1  - ka)/kax(l-k2)  . 

This  expression  is  negative  for  0 < k < Y”1^2*  positive  for  y"1^2  < k 
and  negative  for  k > 1. 

We  find  from  differentiating  Eq.  86  that 


[vr]  - -CY-U^/YP^jU-k1)  , 

which  is  negative  for  k < 1 and  positive  for  k > 1. 
These  results  are  shown  in  Fig.  7 with  k < 


con- 
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(104) 

(105 

(106) 

(107) 

(108) 
: 1 
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When  $2  ■ 0 with  a jump  in  ^ at  rc»  no  jumps  occur  in  these  first 
order  spatial  derivatives  and  we  examine  the  second  order  derivatives. 

In  this  case  we  obtain  similar  forms,  with 

[urrj  - (Y-Dq*  /aJO-ka)  , (110) 

2 

IP„]  - <Y-»L/rl•1(X■*,)  . (HD 

I»„l  “ - 4r_CY',-k!)/k»1Cl-kI)  , (111) 

and 

Ivrr]  - - (Y-l)^  /Ylcp^jCl-k2)  , (113) 

2 

whose  behaviors  are  easily  interpreted. 

These  results  are  summarized  in  Fig.  8 , again  with  k < y-1/2. 

i 

c)  Shock  Wave  Surface 

When  k > 1 with  the  interface  being  a shock,  the  variables  obviously 
jump  across  the  shock,  so  the  following  shock  equations  must  be  satisfied 


a2  " ai C1  + m2)  CyMa  * ^T3  (yTT)2/Mi  • 

P m n M2  - Xli.) 

P2  P1  (y+1  Y+iJ 

“a  ' ui  * Yir*i(M  " M"lj  ’ 


1 (¥  * «-) 


(114) 

(115) 

(116) 
(117) 


b'ssfe  iaSaii.  l ^ 


--tummi . 


where 


M - (i^-up/aj  - k 


(118) 


Fig.  8.  Property  variations  as  a function  of  radius  when  both  q and  qr 

are  discontinuous  at  the  kernel  edge  and  that  edge  moves  as  a 

-1/2 

subsonic  wave  across  particle  paths  with  < < Y . (a)  uy  vs  r 

(b)  p„  vs  r,  (c)  o vs  r and  (d)  v,  vs  r. 
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When  Eqs.  115-118  are  differentiated  following  the  shock,  we  obtain 


and 


where 


^2  . foc_  M*  . I4l£l  ♦ II-  p M 42.  , 

ot  Iy+1  y+1 J 6t  Y+l  J1  it  * 


(119) 


6u,  6u.  , 6a.  .u 

W * St1  + tFT  cm_m‘1)  It1  + C1+M“2)  al  St  * 


(120) 


6v, 


'2  2 

W-mrZ 


+ M‘2)^-  ^iM‘3§ 


(121) 


$r"  (*c  “ sr)/ai " f'c-ui)ar  w-  * 


(122) 


6p 


ST  - -YP^u  + juj/r^  + Crc-Ul)Pj 


6u 


StT  " "Vr,  + <Vul)ui 


(123) 

(124) 


6v 


ST*  vl(ur1  -W  * t'cY’r,  ’ 


(125) 


6a 


sr-  - Ilri*i‘ur1  ♦ W * <vui,*r1 


(126) 


«P2 


ST  • * Ju2/rc>  * tY-l)q2/Vj  ♦ (ic-«2)prj  . (127) 


6u. 


St Vr2  + Crc“u2)ur2  • 


(128) 


ST  " v2(ur2  + *u2/rc>  ♦ Crc-u2)V|.2 


(129) 


*k 7'“' 


’ " * ■-■ 
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Solving  Eqs.  127  and  128  for  u and  p , we  obtain 

r2  r2 


/ 6u  6p,  v 

V * {CVU23  It  + V2  It  “ CY“1)q2  + a2jVrc//{(VU2)4"a2}  Cl30) 


and 


pr2  * )(VU23 


«p 

IT  “ CY‘1)q2/v2  + YP2ju2/rc 


^2  srt/icvu2),-*2 


From  Eqs.  129  and  86,  we  obtain 


6v, 


I ST*  - Vur2  * Ju2/rc> 


}/ftc-u; 


} 


(131) 


(132) 


and 


(133) 


It  should  be  noted  that  Eqs.  130-133  can  be  expressed  in  terms  of 
(1)  q2,  (2)  rc  and  its  time  derivatives  and  (3)  the  flow  variables  and  their 
spatial  derivatives  evaluated  just  ahead  of  the  shock.  However,  the  re- 
sulting expressions  are  quite  long,  so  have  not  been  included  here. 

The  explicit  method  of  integration  in  the  CLOUD  program  does  not 
properly  handle  the  jumps  and  corners  in  the  variables  which  result  from 
discontinuities  in  q or  ^ at  the  kernel  edge.  As  a consequence,  the  com- 
putation becomes  unstable  and  soon  "blows  up."  However,  if  q is  sufficiently 
smooth,  with  both  4 and  ^ equal  to  zero  at  the  kernel  edge,  these  jumps  and 
corners  disappear.  Then  CLOUD  seems  to  work  very  nicely  so  long  as  reaction 
may  be  neglected. 
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7.11  Status  at  the  End  of  Grant  Period 

Calculations  for  several  time  lines  beyond  the  initial  time  line  from 
the  CLOUD  program  have  been  made,  assuming  in  the  implicit  program  that  in- 
tegration along  the  particle  direction  as  well  as  along  the  K+  and  K 
directions  can  be  adequately  represented  by  the  arithmetic  average  of  values 
for  the  two  end  points.  Examination  of  the  computed  values  for  uy  and  pr  near 
the  origin  indicated  that  these  variables  did  not  seem  to  be  behaving  pro- 
perly. The  reason  appears  to  be  that  for  the  K+  and  K-  directions  the  pre- 
vious integration  assumption  is  inadequate,  as  we  will  now  show. 

Let  us  consider  the  integral  of  * u-K  a near  the  origin.  The 
behavior  of  u vs  r for  time  lines  t and  tn+1  are  approximately  as  shown  in 
Fig.  9.  If  we  use  average  values  corresponding  to  the  end  points,  we  get 

rCl.nU)  - rC2.nl  - jMM+l^MJ. 

. K-  I»(l,n»l)«C2,n)J  j[tn<1-tn]  . (134) 

However.  uCl.n.l)  - u(l,n)  « 0,  so  the  contribution  from  u is  “'^'"'-(t^j-t,,) . 
which  is  simply  the  average  value  of  u multiplied  by  the  time  increment. 

Next,  we  consider  the  dashed  line  which  approximates  the  path.  Clearly,  the 
integral  of  u along  this  path  is  greater  than  u C?.fffl(tntl-tn) , so  we  require 
a better  approximation  to  u. 

We  obtain  a better  approximation  if  we  assume  u may  be  represented  by 
a function  containing  two  groups  of  terms.  The  first  group  represents 
how  u varies  along  the  t ■ tnline;  the  second  group  represents  how  much  u 
has  changed  in  time  (t  - tn)  from  its  value  at  tR.  Thus  we  write 


VELOCITY 


0.3 


0.2 


0.0 


RADIUS  , r ( cm  x I0"3 ) 

Fig.  9.  Behavior  of  particle  velocity  versus  radius  near  the  origin  for 
tines  c_  and  t . . . 


.so- 


il - uci.n)  + urCl,n)lr-r(l,n)J  + urr(l.n)^r— 

♦ j«£l»n)  ♦ ur(l»n)[r-rCl,n)3  + ufr(l,n)  |tt-tn3  , (13S) 

where  u ■ jj-  end  r-r(l,n)  is  given  to  first  order  in  t by 

r-r(l.n)  - -feC^,T  lO>SU  [t-t  ] . (136) 

n+1  n 


Expressions  similar  to  u can  be  written  for  a or  for  other  variables  when 
6K‘  fiK- 

integratin.  and  . A similar  concept  can  be  applied  when  in- 
tegrating in  the  K*  direction.  Changes  in  the  computer  program  to  incor- 
porate these  concepts  were  underway  at  the  end  of  the  grant  period. 

A second  conversion  also  was  underway,  for  the  University  computing 
facility  was  in  transition  from  an  IBM  360  system  to  a CDC  system.  Tho 
switch  to  a different  machine  requires  conversions  of  programs  written 
specifically  for  IBM  machines. 

As  this  account  has  tried  to  indicate,  difficulties  have  arisen 
during  the  period  of  the  grant.  Generally  we  have  been  able  to  identify  the 
cause,  to  analyze  why  it  has  given  us  difficulty  and  to  devise  a method  to 
overcome  each  difficulty.  In  conclusion,  while  we  have  not  yet  developed  a 
computer  program  which  successfully  used  an  implicit  technique  to  solve  our 
problem,  we  still  believe  our  approach  can  succeed  and  intend  to  pursue  the 
matter  further. 
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VIII  DIRECT  INITIATION  OF  DETONATION  BY  NON-IDEAL  BLAST  WAVES 

This  work  is  fully  described  in  AFOSR  TR -77-084 2,  Therefore,  only  a 
brief  description  of  the  work  will  be  given  here.  The  reader  is  referred 
to  the  full  report  for  more  detailed  description. 

The  purpose  of  this  study  was  to  examine  the  process  of  direct  ini- 
tiation for  a non-ideal  blast  wave  source  for  reasonably  realistic  kinetics 
in  the  gas  phase  region  surrounding  the  source.  The  kinetics  was  modeled  as 
Arrhenius  kinetics  for  the  hydrogen-oxygen  reaction  but  it  was  assumed  that 
the  ignition  delay  time  was  infinite  outside  the  temperature  range  of  1000  K 
to  2700  K.  This  is  a reasonably  realistic  model  of  the  actual  kinetics  be- 
cause below  1000  K the  delay  to  explosion  increases  very  rapidly  with  drop- 
ping temperature  while  above  about  2700  K the  combustion  reaction  in  the 

i 

hydrogen -oxygen  system  becomes  endothermic  because  of  the  excessive  amount 
of  dissociation  produced  by  the  high  temperature.  Thus,  in  the  hydrogen-oxy- 
gen system  one  would  expect  to  see  an  explosion  induced  in  the  high  tempera- 
ture gas  only  over  the  range  of  about  1000  to  2700  K. 

The  calculations  were  performed  as  follows.  Existing  CLOUD  output 
for  various  cases  of  bursting  spheres  and  ramp  addition  of  energy  plus  a 
number  of  new  calculations  for  specific  cases  were  used  to  describe  the  blast 
wave  flow  from  a typical  non-ideal  source  (either  a bursting  sphere  or  a 
region  where  the  energy  is  added  as  a ramp  addition  of  energy).  This  last 
addition  quite  closely  approximates  spark  addition  of  energy.  These  back- 
ground flows  were  all  calculated  in  dimensionless  coordinates.  Thus  multi- 
plication by  the  proper  dimensional  quantities  allows  one  to  change  the  actual 
size  of  the  source  region  without  recalculating  the  entire  flow  field.  A 
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program  was  constructed  which  examined  the  time,  temperature  and  pressure 
history  of  cells  number  52  through  130  for  these  different  runs.  In  each 
case,  along  each  Lagrangian  cell,  which  represents  a particle  path,  hydro- 
gen-oxygen  kinetics  was  assumed  and  the  accumulated  delay  to  explosion 
was  calculated  using  an  experimentally  determined  Arrhenius  equation  for 
the  delay.  In  one  or  more  of  these  cells  the  accumulated  delay  to  explosion 
had  a minimum  time  and  this  was  assumed  to  be  the  point  in  the  spherical  cell 
surrounding  the  Blast  source  which  first  exploded.  The  calculated  dimension- 
less time  to  explosion  was  converted  to  a real  time  fractional  delay  to 
explosion  by  using  the  size  of  the  blast  source  as  a scaling  par.  . Mr.  As 
the  source  gets  smaller  for  any  fixed  energy  density  source  the  total  energy 
in  the  source  gets  smaller  and  at  the  same  time  the  real  accumulated  frac- 
tional explosion  delay  time  gets  smaller  eventually  reaching  a value  of  one. 
Under  these  circumstances  a smaller  sized  source  will  not  produce  direct 
initiation  of  detonation.  Thus,  one  can  systematically  study  the  effect  of 
energy  addition  rates  and  of  energy  density  as  well  as  total  energy  of  the 
source  region  on  the  initiation  behavior  for  an  hydrogen-oxygen  detonation 
system.  This  was  done  and  the  results  show  that  for  sources  with  low  energy 
density  the  source  energy  required  to  produce  detonation  rises  to  infinity 
while  for  high  energy  density  there  is  an  asymptotic  minimum  value  of  source 
energy  which  is  required  for  initiation.  This  behavior  agrees  quite  well 
with  the  recent  description  of  initiation  that  has  been  presented  by  Professor 
John  Lee  of  .McGill  University.  Bven  though  this  study  showed  ignition  delay 
times  and  initiation  energies  which  were  well  below  what  one  finds  experi- 
mentally for  this  type  of  source  behavior  the  trends  are  qualitatively  cor- 
rect. The  large  quantitative  difference  (our  energies  are  approximately  a 
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factor  of  100  lower  than  the  experimental  energies)  can  be  related  to  the 
fact  that  direct  initiation  only  occurs  when  there  is  sufficient  time  to 
produce  transverse  waves  of  sufficient  amplitude  to  cause  self- sustenance  of 
the  wave.  The  simple  theory  in  the  report > however,  predicts  the  time  at 
which,  direct  homogeneous  initiation  would  occur  and  has  nothing  to  say  about 
the  development  of  transverse  wave  structure. 

DC  OTHER  WORK.  OF  INTEREST 

There  has  been  a considerable  amount  of  effort  expended  over  the 
period  of  the  grant  to  study  the  effect  of  source  behavior  on  the  nature  of 
the  blast  wave  produced  when  the  source  is  a low  energy  density  source  simi- 
lar to  that  which  occurs  in  a vapor  cloud  explosion  or  in  the  explosion  of  an 
FAE  device.  Briefly,  systematic  studies  have  been  performed  for  blast  waves 
produced  by  bursting  spheres,  by  the  ramp  addition  of  energy  (which  models  a 
spark)  and  by  constant  velocity  and  accelerating  flames  from  very  low  Mach 
numbers  up  to  Mach  numbers  above  the  Chapman- Jouguet  Mach  number.  As  a re- 
sult of  these  studies  the  behavior  of  the  blast  wave  from  bursting  spheres 
and  from  constant  velocity  flames  can  now  be  estimated  quite  accurately  if 
either  the  initial  conditions  of  the  bursting  sphere  are  known  or  if  the  com- 
bustion system  and  the  flame  velocity  are  known  for  the  constant  velocity 
flame  case.  Additionally,  it  has  been  shown  that  a spherical  accelerating 
flame  will  not  produce  overpressures  in  excess  of  that  produced  by  the 
highest  flame  velocity  in  the  system  when  the  explosion  consists  of  a flame 
traveling  at  that  high  velocity.  This  means  that  calculations  for  constant 
velocity  flames  are  adequate  to  estimate  blast  damage  if  the  maximum  flame 
velocity  is  known.  This  work  has  lead  to  one  M.S.  and  two  Ph.D.  theses  and 
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to  four  papers,  two  of  which  have  been  published  and  two  of  which  have  been 
submitted  to  journals  for  publication  (15,18-23).  Additionally,  two  M.S. 
theses  are  currently  in  the  process  of  preparation.  These  will  cover  the 
experimental  work  on  initiation  delay  times  for  propylene  oxide-nitrogen- 
oxygen  mixtures  and  the  effect  of  inhibitors  and  promoters  on  ignition  de- 
lay times  for  propane-air  mixtures. 
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